use "${clean}newspaper_analysis.dta", clear
keep if north == 1 | south == 1

// recoding Alexandria, VA to the South (technically part of DC until 1844)
replace south = 1 if city == "Alexandria"
replace north = 0 if city == "Alexandria"
replace state = "Virginia" if city == "Alexandria"

// dropping DC as a border area
drop if state == "District Of Columbia"

// standardizing variables
egen county_surban = std(county_urban_percent)
egen county_srugged = std(county_ruggedness)
egen county_spostoffice = std(county_postoffice_area)

// keeping the sample consistent
cap drop sample
gen sample = 1
foreach var of varlist $county_npcontrols {
	replace sample = . if `var' == .
}

// creating some variables
local c = 1865

gen period = year > `c'
replace year = year - `c'

encode state, gen(state_id)

	
// matrices to store the results for graphing
matrix define R = J(5, 3, 0)


* regressions
* ------------

// City FEs
xtset city_id
qui xtreg singular c.year##i.period##i.north ${county_npcontrols} if sample == 1, fe cluster(city_id)
	
qui lincom 1.period // intercept change at cutpoint
matrix R[1,1] = `r(estimate)'
matrix R[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[1,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix R[4,1] = `r(estimate)'
matrix R[4,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[4,3] = `r(estimate)' - 1.96*`r(se)'
	

// State FEs
xtset state_id
qui xtreg singular c.year##i.period##i.north ${county_npcontrols} if sample == 1, fe cluster(state_id)

qui lincom 1.period // intercept change at cutpoint
matrix R[2,1] = `r(estimate)'
matrix R[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[2,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix R[5,1] = `r(estimate)'
matrix R[5,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[5,3] = `r(estimate)' - 1.96*`r(se)'


* graphing 
* --------
svmat R

// temp variable to index estimates
cap drop temp*
gen tempn = _n if R1 != .

gen tempslope = tempn > 3

gen tempfe = mod(tempn, 3)
recode tempfe (1=0) (2=1) (0=.)

 
twoway ///
	(scatter R1 tempn if tempslope == 0 & tempfe != .,	msize(large) mlwidth(none) mcolor("${red}") 	msym(D) yaxis(1 2) ) ///
 	(rcap R2 R3 tempn if tempslope == 0 & tempfe == 0,   lwidth(thick) lcolor("${red}") 	lpattern(dash) yaxis(1 2)) ///
 	(rcap R2 R3 tempn if tempslope == 0 & tempfe == 1,   lwidth(thick) lcolor("${red}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(1 "City FEs" 2 "State FEs") ///
	xtitle("") 	///
	ylab(none, axis(2)) ///
	ylab(-.25(.125).25, axis(1) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(0.5 2.5)) ///
	subtitle("Intercept Shift: 1865 - 1866", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	fxsize(100) fysize(100) ///
	name(intercepts, replace) nodraw
	
	
twoway ///
	(scatter R1 tempn if tempslope == 1 & tempfe != .,	msize(large) mlwidth(none) mcolor("${red}")	msym(S) yaxis(1 2)) ///
	(rcap R2 R3 tempn if tempslope == 1 & tempfe == 0,   lwidth(thick) lcolor("${red}") 	lpattern(dash) yaxis(1 2)) ///
	(rcap R2 R3 tempn if tempslope == 1 & tempfe == 1,   lwidth(thick) lcolor("${red}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(4 "City FEs" 5 "State FEs") ///
	xtitle("") 	///
	ylab(none, axis(1)) ///
	ylab(-0.02(0.01)0.02, axis(2) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(3.5 5.5)) ///
	subtitle("Slope Change after 1865", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	fxsize(100) fysize(100) ///
	name(slopes, replace)  nodraw
	
graph combine intercepts slopes,  scale(1.5) xsize(100) ysize(50)
graph export "${output}FigD5_np_south_1865_altfes.pdf", as(pdf) replace
